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Formation and Collapse of Nonaxisymmetric Protostellar Cores in 

Planar Magnetic Molecular Clouds 

Shantanu Basu^ and Glenn E. Ciolek^ 
ABSTRACT 

We extend our earlier work on ambipolar diffusion induced formation of pro- 
tostellar cores in isothermal sheet-like magnetic interstellar clouds, by studying 
nonaxisymmetric collapse for the physically interesting regime of magnetically 
critical and supercritical model clouds (/io > 1, where /xq is the initial mass- 
to- magnetic flux ratio in units of the critical value for gravitational collapse). 
Cores that form in model simulations are effectively triaxial, with shapes that 
are typically closer to being oblate, rather than prolate. Infall velocities in the 
critical model (/iq = 1) are subsonic; in contrast, a supercritical model (/iq = 2) 
has extended supersonic infall that may be excluded by observations. For the 
magnetically critical model, ambipolar diffusion forms cores that are supercrit- 
ical (yUcorc > 1) and embedded within subcritical envelopes (/ienv < !)• Cores 
in our models have density profiles that eventually merge into a near-uniform 
background, which is suggestive of observed properties of cloud cores. 

Subject headings: diffusion — ISM: clouds — ISM: kinematics and dynamics — 
ISM: magnetic fields — MHD — stars: formation 



1. Introduction 



Magnetic fields play an important role in star formation, especially in the early stages 
of core formation and collapse; measured mass-to-flux ratios of molecular clouds yield an 
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average that is ~ 1 — 2 times greater than the critical value for collapse (Crutcher 1999). 

However, observational biases tend to push toward higher values of measured mass-to-fiux 
ratio (Crutcher 2003, private communication), so that moderately subcritical cloud regions 
arc not ruled out. Dense cores within molecular clouds arc the sites of star formation, with 
detected infall up to ~ 0.5 Cg ~ 0.1 km s~^ on scales < 0.1 pc, (e.g., in L1544, Tafalla et al. 
1998; Williams et al. 1999), where Cg is the isothermal sound speed. Ciolek & Basu (2000) 
have fit the main features of the observed velocity and density profiles in L1544, modeling 
it as an axisymmetric supercritical core embedded in a moderately subcritical envelope. 

Axisymmctry is clearly an idealization to real cores. Observations suggest a typical 
projected axis ratio of 0.5 for cores (Myers et al. 1991), and deprojections of the distribution 
of shapes imply intrinsically triaxial but nearly oblate cores (Jones, Basu, & Dubinski 2001). 
Polarized emission measurements from dense cores also imply triaxiality (Basu 2000). More 
generally, detailed submillimeter maps of star-forming regions reveal significant irregular 
structure and multiple cores (Motte, Andre, & Neri 1998). Theoretical nonaxisymmetric 
magnetic models of the collapse and fragmentation of a single core were presented by Naka- 
mura & Hanawa (1997) and Nakamura & Li (2002), without and with ambipolar diffusion, 
respectively, using the magnetic thin-disk approximation (Ciolek & Mouschovias 1993, here- 
after CM93). The early stages of core formation in a nonaxisymmetric infinitesimally thin 
subcritical sheet (including the effects of magnetic tension but ignoring magnetic pressure) 
were studied by Indebetouw & Zweibel (2001). Here, we also study a planar cloud that is 
perpendicular to the mean magnetic field, focusing on the case of either exactly critical or 
decidedly supercritical cores; these cases are shown to lead to observationally distinguish- 
able outcomes. We again use the thin-disk approximation, which allows for finite thickness 
effects, and explicitly includes the effects of both magnetic pressure and tension. 

2. Physical and Numerical Model 

The fundamental equations we use to model molecular clouds as self-gravitating, partially- 
ionized, isothermal thin sheets or disks have been presented in axisymmetric form in CM93 
and Basu & Mouschovias (1994, hereafter BM94). In this study, the condition of axisymmc- 
try is no longer employed, and instead we model clouds as thin sheets of infinite extent [with 
vertical half-thickness Z{x, y, t)] in a cartesian coordinate system {x, y). Magnetohydrostatic 
equilibrium in the z-direction (i.e., balance of thermal- pressure, gravitational, and magnetic 
pinching forces) is assumed at all times. The evolution is followed in a square region using 
periodic boundary conditions. To find the gravitational field components in the sheet, we 
calculate the gravitational potential '>p{x,y,t). For this situation, solving Poisson's equation 
(see eq. [29] of CM93) relates to the column density aa(x,y,t) (— J^^ pa(x,y,t)dz; is 
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the mass density) by 

= -2T,Gr[a^]/K, (1) 

where T[f \ is the Fourier transform of a function /, kz {— [kl + kyY^'^) is the wavenumber, 
and G is the gravitational constant. In our governing equations, the effects of magnetic 
pressure and magnetic tension are both included. The magnetic potential '^{x^y^t) that is 
used to determine the x- and y- components of the magnetic field at the surface of the sheet 
(necessary to calculate magnetic tension forces) is also given by equation (1), by substituting 
%l) with ^, and 27rG'o"n with — (i?2,cq — -Brcf); for details see CM93. Bz^cq{x^ y, t) is the vertical 
magnetic field strength in the sheet, and i?ref is the uniform magnetic field of the background 
(reference) state, as well as the assumed constant magnetic field as \z\ —>■ oo. In the present 
study, we neglect magnetic braking due to a finite density external medium (BM94), and 
dust grains and UV ionization (CM93, Ciolek & Mouschovias 1995). 

In our formulation, velocities are normalized to the isothermal speed of sound Cg = 
0.188 (T/ 10 K)^/^ km s~^ (where T is the temperature and we have used a mean molecular 
mass m — 2.33mH, in which mn is the hydrogen atom mass), column densities are in units of 
(7n,o, the uniform column density of the background state. The time unit is to — Cs/27rG'(7n,o, 
the length unit is Lq = c^/2TiGa^fl^ and the magnetic field unit is Bq = 27rG^''^(Tn,o- As 
discussed in CM93 and BM94, model clouds arc distinguished by three basic dimensionless 
parameters, namely, the initial mass-to-magnetic flux ratio (in units of the critical value for 
gravitational collapse) /xq = Bq/B^cI = 27rG'^/^an,o/-Brcf, the initial dimensionless neutral- 
ion collision time fni.o = Tni/to; and the normalized external pressure acting on the disk, 
Pext = -Pext/ (f G^cTn.o)- The iou number density rii is calculated using the scaling law rii oc rh/"^ 
(as done by BM94; nn is the neutral number density), which is a reasonable approximation 
for Un < 10^ cm~^ (e.g., Ciolek & Mouschovias 1998). For a background number density 
JT-n.o = 3 X 10^ cm~^ and Pcxt = 0.1, we find that crn,o = 5.98 x 10~^ g cm~^ (therefore 
Nn,o = CTnfi/m = 1.54 X 10^1 cm-2), Lq = 4.57 x 10"^ pc, and to = 2.38 x 10^ yr. 

We use the method of lines technique (e.g., Schiesser 1991); the system of partial 
differential equations of two-fiuid MHD are converted to a system of ordinary differential 
equations (ODE's) through second-order spatial finite differencing, and use of the van Leer 
(1977) advection scheme. Time integration of the ODE's is performed using the implicit 
Adams-Bashforth-Moulton method. Details of our basic technique can be found in Morton, 
Mouschovias, & Ciolek (1994). Fast Fourier transform (FFT) routines are used to solve for 
the gravitational and magnetic potentials obtained through equation (1). The computational 
region has extent 4AT,m on each side, where AT,m = 47rLo is the wavelength with maximum 
growth rate of gravitational instability in a non-magnetic (thermal) infinitesimally thin disk. 
Models presented in this paper are run on a uniform grid of 128 x 128 points. 
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3. Results 

We present the evolution of two model clouds; one model is exactly critical, with fiQ = 1, 
and the other is supercritical with = 2. Both models have fni,o = 0.16 (as adopted in the 
standard model of BM94) and Pext = 0.1. Since the background (reference) state is charac- 
terized by a uniform column density (Jn,o and magnetic field Bj-^iZ (where Bj-^i — Bq/ hq) the 
gravitational and magnetic forces in the reference state are identically zero. To initiate evo- 
lution, we superimpose a set of random, small-amplitude (the rms is 3% of the background) 
column density pcrtTirbations Sa^- The magnetic field is also perturbed, with SBz^eq ^cTn, 
to maintain flux-freezing in the initial state of each model. The spectrum of perturbations 
is flat (white noise) , with damping so that wavelengths of twice the grid spacing and smaller 
are negligible; we choose this particular type of spectrum so as not to preferentially excite 
modes with the maximum growth rate. The evolution of each model is followed until the 
maximum column density 10 (Tn,o- Beyond this enhancement, gravitational instability can- 
not be spatially resolved, and the evolution is relatively very rapid in a small region near 
the density peaks, making a simulation of the larger cloud impractical even at higher resolu- 
tion. Since the vertical balance along field lines is primarily between gas pressure pn Cg and 
self-gravitational pressure f Gcr^, the density has increased by a factor ^ 100. 



3.1. Critical model (/xq = 1) 

Figure 1 exhibits color images of the normalized column density (Tn{x,y)/(Tn,Q and the 
local mass-to- fiux ratio fJ>{x,y) — 2nG^^^an/B2^e(i over the entire computational region for 
the critical model at its final output time, t — 133.9^0- The time is rather long due to the 
very small-amplitude perturbations and because magnetic forces exactly balance gravity in 
the flux-freezing limit, requiring instability to be initiated by random ion-neutral drifts. Two 
major density peaks have formed by this time. The column density structure shows that the 
lowest density contours are the most irregular and elongated. Henceforth, we define a dense 
core to be the region within the an/a^fi = 2 contour, corresponding to the observationally 
significant threshold rin ~ 10^ cm~^. The masses of the two dense cores in this simulation 
are 5.3Mq (lower left of Fig. la) and 4.2Mq (upper right of Fig. la), respectively, using 
the standard values of units given in § 2. For comparison, the enclosed mass in our periodic 
computational domain is 151Mq. The shape of the cores are generally nonaxisymmetric, a 
result that was first seen in the magnetically subcritical planar ambipolar diffusion models 
of Indebetouw & Zweibel (2000). The mean value of Z = a^/{2p^) implies one oblate core 
with relative axis lengths 0.15:1:1, and another that is triaxial with relative lengths 0.1:0.7:1, 
making it more nearly oblate than prolate. The cores formed in this model always have their 
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shortest axis aligned normal to the sheet (parallel to the mean magnetic field) . Although one 
core is effectively oblate in this particular model realization, our simulations reveal they are 
most commonly triaxial and more nearly oblate than prolate, in accordance with statistical 
analyses of observed cores (Jones et al. 2001). 

Velocity vectors of the neutrals are also displayed in Figure la. Calculation of the 
rms speeds of neutrals (vn,rms) and ions (fi,rms) in the entire computational region yield 
■i^n.rms = 0.114 Cg = 1 .14v[^„^^s] iohs lag behind neutrals because the motions are ultimately 
gravitationally driven. Furthermore, all of the infall speeds in this model are subsonic. The 
maximum infall speed is |t'n| = 0.54 Cs, and speeds of this order occur only within the cores. 
As we shall see in § 3.2, subsonic infall is a distinct feature of ambipolar-diffusion initiated 
collapse in clouds with /xq < 1- This trend is consistent with multi-line studies that find 
subsonic gas speeds in dense protostellar cores (Tafalla et al. 1998; WiUiams et al. 1999), 
and with axisymmetric models of core formation in subcritical clouds (Ciolek & Basu 2000). 
Figure lb shows another interesting result: the cores are supercritical (//core > o-i^d are 
surrounded by magnetically subcritical envelopes (fJ-env < This is a natural consequence 
of ambipolar diffusion, which redistributes mass in magnetic flux tubes (Mouschovias 1978). 
Due to the initial precise balance between gravitational and magnetic forces in this model, 
evolution occurs as ambipolar diffusion effects a drift of mass through essentially stationary 
magnetic field lines. In time, this leads to the formation of both supercritical cores and a 
subcritical envelope. This can also occur in clouds with /^o slightly above unity. 

Figure 2 shows the column density c^/a^fi, vertical magnetic field B^cqZ-Bo, and x- 
component of neutral velocity Vn,x/ Cs along the x-axis for a line that cuts through one of the 
cores shown in Figure 1. The core stands out as a well-defined high density region with an 
eventual merger into a lower-density background (the remnant of the initial uniform state) 
that surrounds it; this is reminiscent of the mid-infrared absorption maps of dense cores 
made by Bacmann et al. (2000). The magnitude of v^^^ increases inward toward the core 
center (but remains subsonic throughout) before dropping to zero at the core center. 

3.2. Supercritical model (/^o = 2) 

Figure 3a exhibits a color image and contour plot of cr^/cr^fi, as well as velocity vectors 
for the supercritical model at its final output time, t — 17.6 to- The time is significantly 
lesser than in the previous model due to the relative ease of gravitational instability in a 
supercritical cloud. Because the evolution is more rapid and dynamical in the supercritical 
model, ambipolar diffusion doesn't have as much time to operate. Hence, the mass-to-flux 
ratios of the high-density cores in this model (/^core — 2.1) are only slightly greater than the 
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Fig. 1. — Critical model (/io = 1). The data are shown when the maximum column density 
~ lOdn^o- (fl; Left.) Image and contour plot of normalized column density crn(x, ?/)/(Tn,o- 
The contour lines are spaced in multiplicative increments of 2^/^. Also shown are velocity 
vectors of the neutrals; the distance between tips of vectors corresponds to a speed 0.5 Cg. {b, 
Right.) Image and contours of fi{x, y), the mass-to-fiux ratio in units of the critical value for 
collapse. Regions with fj, < 1 are black. The contour lines are spaced in additive increments 
of 0.1. 
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initial value. The core shapes are again nonaxisymmetric, either near-oblate or decidedly 
triaxial. For this model we find i'n,rms — 0.222 Cg = 1.05 T;i^i.ms; these results differ from the 
critical model due to the more dynamical evolution and greater ability of neutrals to drag 
ions and magnetic field with them. The maximum infall speed is 1.19 Cg, and speeds of 
this order typically occur within distances ^ 0.1 pc from the core centers. In contrast to 
the critical model, the infall speeds within the cores in this model are often supersonic, and 
significant motions 0.5 Cg are also seen throughout the cloud, and well outside the cores. 
Figure 36 shows Cj^/a^fi, B^^eq/Bo, and Vn,x/cs along the x-axis for a line passing through 
one of the cores in Figure 3a. The collapsing core again has a density profile that eventually 
merges into the background. 

Our model shows that the extent of the supersonic fiow region is well within the res- 
olution capability of current observations (scales < 0.1 pc, n-a ~ 10^ — 10^ cm"^). Since 
supersonic infall has not been detected over these length scales and densities, and models 
with /iq > 2 would have even greater infall speeds, this suggests that molecular clouds that 
are supercritical by a factor > 2 are incompatible with observations of protostellar cores. 

4. Conclusions 

Ambipolar diffusion leads to a nonuniform distribution of mass-to-fiux ratio, in a nat- 
ural extension of the process described by Mouschovias (1978). Stars form preferentially in 
the most supercritical regions. Core shapes are somewhat triaxial, and usually more nearly 
oblate than prolate. The core column density eventually merges into a near-uniform back- 
ground value. In the critical (/io = 1) model, a surrounding region is established which is 
mildly magnetically subcritical (due to flux redistribution) and infall motions both inside 
and outside cores are subsonic. Conversely, the supercritical (/^o = 2) model exhibits super- 
sonic motions within cores, and extended rapid motions outside them. The critical model 
requires a significantly longer time to develop gravitational instability; however, we caution 
that the growth time for both models are likely upper limits due to the possibility of nonlin- 
ear perturbations in more realistic situations. We also note that all motions in these models 
are fundamentally gravitationally driven; the neutral speeds are typically greater than those 
of the ions. 

We thank Craig Markwardt for providing some key IDL routines. SB is supported by 
a grant from the Natural Sciences and Engineering Research Council of Canada. GC is 
supported by NASA grant NAG 5-7589 to the New York Origins of Life Center (NSCORT). 
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Fig. 2. — Physical quantities in the //q = 1 niodel, along a line parallel to the x-axis that cuts 
through a supercritical core (shown in Figure 1) centered a,t x — 0.64 AT,m, y — 0.70 AT,m- 
Solid curve: column density (Tn/crn,o- Dashed curve: magnetic field strength S;j_eq/-Bo- Dotted 
curve: x-component of the neutral velocity, Vn,x/cs- 



- 10 - 




-2 -1 1 2 0.8 1.0 1.2 1.4 1.6 1.8 2.0 



Fig. 3. — Supercritical model (/iq = 2). The data are shown when the maximum column den- 
sity ~ lOcJn^o- (fl) Left.) Image and contour plot of normalized column density (Tn(x, |/)/(Jn,o- 
The contour lines are spaced in multiplicative increments of 2^/^. Also shown are velocity 
vectors of the neutrals; they are normalized in the same way as in Fig. la. {b, Right.) 
Physical quantities along a line parallel to the x-axis that cuts through a supercritical core 
centered at x = 1.42 AT.m,?/ = — 0.33AT,m- Solid curve: column density cTn/cTn.o- Dashed 
curve: magnetic field strength Bz^cq/ Bq- Dotted curve: x-component of the neutral velocity, 



